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ABSTRACT 

This  thesis  describes  the  development  of  a  Fortran  IV 
computer  program  for  finding  shock-induced  stresses  in  a 
three-dimensional  piping  system.   Discretized  equations  of 
motion  are  formed  by  the  finite  element  method.   Input 
motions  are  support  translations  specified  by  response 
shock  spectra.   Maximum  responses  in  individual  modes  and 
resulting  octahedral  shearing  stresses  at  selected  points 
are  found  for  shock  motions  in  three  orthogonal  directions, 
Extreme  stresses,  based  on  the  square  root  of  the  sum  of 
the  squares  of  the  modal  stresses,  are  estimated  for  each 
input  direction.   An  example  problem  is  analyzed  to  demon- 
strate the  use  of  the  program. 
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I.   INTRODUCTION 

With  the  greater  demand  for  electrical  power,  and  the 
diminishing  fossil  fuel  supply,  there  has  been  an  increase 
in  the  construction  and  use  of  nuclear  power  generating 
plants.   Correspondingly,  there  has  been  greater  concern  for 
the  safety  of  the  power  plant  during  an  earthquake. 

Also,  the  reduction  in  the  size  of  the  Navy  has  led  to 
an  increased  emphasis  on  the  integrity  of  the  internal 
systems  of  a  naval  vessel  under  shock  loads.   Thus,  there 
has  been  a  growing  engineering  interest  in  finding  a  rapid 
and  accurate  method  of  analysis  to  determine  the  shock- 
induced  stresses  in  a  complex  continuous  piping  system. 
This  thesis  describes  and  presents  a  computer  program, 
written  in  Fortran  IV  computer  language,  to  find  the 
stresses  in  a  piping  system  responding  to  shock. 

Complex  piping  systems  can  be  analyzed  by  using  dis- 
cretized  models.   In  using  the  discretization  technique,  the 
continuous  piping  system  is  modeled  by  finite  size  elements 
connected  at  nodes.   After  the  system  has  been  subdivided, 
the  elastic  and  inertial  matrices  of  each  element  can  be 
found  and  assembled  to  form  the  elastic  and  inertial 
matrices  of  the  system  [1,2]. 

A  considerable  volume  of  material  is  available  to  assist 
in  analyzing  structures  and  systems  responding  to 
earthquakes  [3,^,5].   Two  of  the  techniques  used  in  earth- 
quake response  analysis  are  modal  analysis  with  a 
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discretized  model  and  transfer  functions  [4,5].   The  shock 
input  for  earthquakes  can  be  specified  by  using  time 
history  [3]>or  shock  spectrum  [4,5]. 

The  Navy  uses  a  technique  called  the  Dynamic  Design- 
Analysis  Method  (DDAM)  to  determine  the  shock  response  of 
shipboard  equipment  and  systems  [6].   The  DDAM  is  a  modal 
analysis  technique  using  shock  spectra  to  specify  the  shock 
inputs. 

The  shock  spectrum  presents,  as  a  function  of  the  system 
natural  frequency,  the  maximum  response  of  a  one  degree-of- 
freedom  system  responding  to  a  shock  motion.   The  shock 
spectrum  is  generally  presented  in  terms  of  response  velocity 
or  acceleration. 

Several  graduates  of  the  Naval  Postgraduate  School  have 
written  theses  on  determining  natural  frequencies  and  mode 
shapes  for  piping  systems.   Fink  [7]  developed  a  program  to 
analyze  planar  piping  systems,  including  out-of-plane 
bending,  using  transfer  matrices.  Baird  [8]  presented  the 
necessary  theory  to  expand  Fink's  work  to  a  three-dimensional 
piping  system.   Rudolf  [9]  developed  a  program  that  deter- 
mines the  frequencies  of  a  general  three-dimensional  piping 
system.   Because  the  transfer  matrix  does  not  yield  a  dis- 
cretized model  of  the  structure,  there  is  no  ready  means 
for  using  the  frequency  and  mode  shape  data  to  find  shock- 
induced  modal  responses.   On  the  other  hand,  the  finite 
element  technique  provides  a  consistent  discretized  model 
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which  may  be  used  for  finding  frequencies,  mode  shapes,  modal 
responses  to  shock  inputs,  and  resulting  stresses.   For  this 
reason  the  finite  element  method  was  chosen  for  this  thesis. 

Shock  spectra  and  mode  participation  factors  are  used 
to  determine  the  modal  responses.   At  selected  points  the 
modal  octahedral  shearing  stresses  are  combined  by  scalar 
addition  of  the  distortion  energy  to  estimate  extreme 
stresses . 
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II.   THEORY 

Material  and  structural  linearity  and  material  local 
homogeneity  and  isotropy  are  assumed  in  developing  the 
mathematical  model. 

A.   MODEL  LIMITATIONS 

There  are  few  inherent  restrictions  on  the  capabilities 
of  the  finite  element  method  to  model  complex  details  of 
piping  configurations.   Considerations  of  the  quantity  of 
input  and  output  data,  program  length,  core  storage  capac- 
ity, computing  time,  and  cost  do  impose  practical  limita- 
tions.  The  limited  time  available  in  developing  this 
thesis  has  necessitated  the  following  restrictions  on  the 
model  used  here. 

The  system  consists  of  straight  lengths  of  constant- 
section  pipe  (elements).   Each  element  centerline  is  parallel 
to  one  of  the  three  mutually  orthogonal  global  refer- 
ence axes,  and  90°  bends  are  replaced  by  fictitious  exten- 
sions of  the  tangent  sections  to  the  intersection  point  of 
the  centerlines.   Effects  of  added  mass  contributed  by 
external  lagging  are  represented,  but  no  provision  is  made 
for  added  mass  due  to  valves  or  fittings.   Pipe  hangers 
furnishing  uniaxial  restraint  in  a  global  direction  may 
also  be  represented.   Bending  action  is  described  by  the 
Euler-Bernoulli  beam  theory,  neglecting  shear  deformation 
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and  rotatory  inertia.   All  ends  of  the  configuration  are 
treated  as  fixed.   The  shock  input  motion  is  a  translation 
of  the  base  along  a  global  axis. 

B.   FINITE  ELEMENT  METHOD 

Once  the  system  has  been  subdivided  into  elements, 
the  elastic  and  inertial  characteristics  are  determined. 
The  element  stiffness  matrix  can  be  determined  by  many 
different  techniques,  whereas  the  inertial  matrix  is  gener- 
ally determined  by  using  shape  functions  and  integrating 
along  the  length,  [1,2]. 

1.   The  Element  Displacement  Vector 

Consider  an  element  with  displacement  components  at 
the  nodes  (ends)  as  shown  in  Fig.  1. 


FIG.  1 

ELEMENT  GENERALIZED  DISPLACEMENTS 


The  clamped  ends  of  the  configuration  are  treated  as 
joined  to  a  rigid  base. 
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Here,  the  double  arrowhead  vectors  are  positive  rotations 
in  accordance  with  the  right-hand  rule.  The  displacement 
vector  is 


T 
q  =  [u1  u2  ...  u12] 


From  Fig.  1,  components  u,  and  u„  are  longitudinal  dis- 
placements, components  u~ ,  u_,  u,-,  Ug,  Ug,  uQ ,  u,,  and  u,« 
are  bending  displacements  and  rotations  in  the  1-2  and  2-3 
planes,  and  components  Ul  and  u..  Q  are  twisting  rotations. 
2.   Element  Stiffness  and  Inertial  Matrices 

Przemieniecki  [2]  forms  the  element  stiffness  ma- 
trix by  solving  the  differential  equations  of  the  element 
displacements.   Appendix  A  reproduces  the  element  stiffness 
matrix  with  Przemieniecki' s  transverse  shear  deflection 
factor  (<{>)  set  equal  to  zero. 

To  determine  the  element  inertia  matrix,  an  expres- 
sion  for  the  element  kinetic  energy  T   is  formed,  rotatory 
inertia  due  to  bending  being  neglected. 


el  -^  .T      T1  itp 

Te    =   j  Ci    f   q      N"    N   q   d£    =   ~  q1    Me   q 

0  " 


where  i   is  the  element  length,  C  is  a  constant,'  Me  is  the 
element  inertia  matrix,  N  is  the  shape  function  row  vector, 
and  %    is  a  dimensionless  length  coordinate. 
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Therefore 


Me  =  Ci    f      NT  N  d£ 


(1) 


Since  the  local  element  axes  coincide  with  the  cross 
section  principal  axes,  the  inertia  matrix  may  be  deter- 
mined from  2x2  and  4x4  element  sub-matrices.   Consider  the 
following  motions  of  the  element: 
a)   longitudinal: 

q^  =  [u]_  u7]T,    N1  =  1  -  5,   N?  =  e 

and  C  =  pA  where  A  is  the  cross-sectional  area  of  the 
element  and  p  is  the  density.   Prom  Eq.  1 


pA£ 


2   1 
1   2 


(la) 


b)   twist: 


qT   =    [U2|    u1Q]    ,       N4    =    (1   -    5)1,      N1Q   =    ?! 

and  C  =  J  where  J  is  the  mass  moment  of  inertia  per  unit 
length   about  the  pipe  axis.   From  Eq.  1 


£-£ 


2   1 
1   2 


(lb) 


c)   bending  in  the  1-2  plane: 


2,o,3 


tel2 


=    [u2   u6   u8   u12]    ,      N2    =    1-35   +2C 


N6    =    £(C-2C2+C3),    Ng    =    3S2-2£3,    N12    =    ft(-€2+C3), 
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and  C  =  pA.   From  Eq.  1 


-B12 


pA£ 
T20 


■13i 

■3*2 


156  22£  5^ 

22£  kl2  13£ 

5^  13  ^  156   -22)1 

-13*  -3S-2  -22£  *U2 


(lc) 


Bending  in  plane  1-3  results  in  the  same  numerical  coeffi- 
cients as  Mm0  .   Because  the  positive  senses  for  the 

— Did 

rotation  angles  are  reversed,  rows  2  and  4  and  columns  2 
and  *J  are  multiplied  by  minus  one.   These  submatrices  are 
assembled  to  form  the  element  inertia  matrix.   This  matrix 
is  displayed  in  Appendix  A. 

C.   ASSEMBLED  EQUATIONS  OP  MOTION 

The  governing  equations  of  motion,  in  matrix  notation, 
for  undamped  free  vibration  [2]  with  a  fixed  base  are 


M  q  +  K  q 


(2) 


where  q  is  the  assembled  displacement  vector  and  M  and  K 

are  the  assembled  inertial  and  stiffness  matrices.   M  and 

e      e 

K  are  referred  to  a  global  coordinate  system  and  M  and  K 

are  referred  to  local  coordinates.   The  assembly  process 
consists  of  realigning  the  local  system  to  the  global 
system  and  building  the  M  and  K  matrices.   Reference  [2] 
gives  sample  assembly  techniques. 
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D.   MODAL  ANALYSIS 

Modal  analysis  is  the  process  of  finding  the  response 
motion  of  a  many  degree-of-freedom  system  by  superposition 
of  the  response  motions  in  the  principal  (or  natural)  modes. 
This  technique  is  advantageous  because  the  motions  in  the 
principal  modes  are  independent  of  one  another  (uncoupled). 
To  use  the  method,  it  is  first  necessary  to  find  the  nat- 
ural frequencies  and  mode  shapes  for  the  free  motions 
governed  by  Eq.  2,  i.e.,  to  solve  the  eigenvalue  problem 

K  v  =  w2  M  v. 

2 

Finding  eigenvalues  (to  )  and  eigenvectors  (v)  is  a 

standard  problem  of  numerical  analysis.   The  first  step  is 
a  transformation  to  a  single-matrix  problem.   Details  of 
the  transformation  used  here,  which  requires  a  Cholesky 
decomposition  of  the  inertia  matrix,  are  given  on  p.  3^9 

of  reference  [1].   Following  this,  Jacobi  rotations  [15] 

2 
are  used  to  find  the  spectral  and  modal  matrices  (ft   and  V) . 

The  modal  matrix  is  normalized  with  respect  to  the  inertia 

T 
matrix,  i.e.,  V  M  V  =  I,  where  I  is  the  nth  order  identity 

matrix. 

The  shock  analysis  is  based  on  using  velocity  shock 
spectra  to  characterize  the  input  motion  and  mode  parti- 
cipation factors  to  determine  the  responses.   Although  both 
shock  spectra  and  mode  participation  factors  have  been 
extensively  used  in  earlier  studies  [6],  the  finite  element 
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discretization  necessitates  a  new  development  of  working 
formulas.   Details  of  this  development  are  given  in 
Appendix  B. 

Each  shock  input  motion  consists  of  base  translation 
parallel  to  a  global  axis.   It  is  shown  in  Appendix  B  that 
the  mode  participation  factor  b   for  mode  r  is 

T 
b   =  v(r)   M  u  (B.5) 

(r) 
where  v    is  the  modal  eigenvector  and  u  is  the  vector  of 

absolute  displacements  corresponding  to  unit  base  displace- 
ment.  The  maximum  relative  displacements  due  to  mode  r 
response  are 


q    =  v(r)  b  V  /uj  (B.6a) 

—max   —     r  r  r 


where  w   is  the  modal  circular  frequency  and  V  is  the 
spectrum  velocity  at  that  frequency. 

Since  responses  to  shock  inputs  in  each  of  the  three 
global  directions  are  separately  determined,  the  vector  u 
and  the  mode  participation  factors  b  must  be  evaluated 
separately  for  each  direction.   The  modal  masses  m   (see 
Appendix  B)  are  likewise  different  for  each  input  direction, 


E.   GENERALIZED  FORCES  AND  STRESSES 

The  subvector  qe  of  element  peak  modal  displacements 
may  be  found  from  the  system  vector,  taking  into  account  the 
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orientation  of  the  local  reference  axes.   From  this  the 
element  generalized  force  vector  f   is  determined.   This 
force  vector  consists  of  the  elastic  and  inertial 
contributions  for  each  element. 


fe  =  (Ke  -  Wr2  Me)  qe 


(3) 


Fig.  2  shows  the  positive  directions  for  the  components 
of  the  generalized  force  vector. 


fe  =  [f  f  ...  f12iT 


FIG.  2 
ELEMENT  GENERALIZED  FORCES 
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The  stresses  at  the  nodal  points  of  the  element  can  be 
determined  from  the  generalized  forces  by  standard  methods 
of  solid  mechanics  [10]. 
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III.   PROGRAM  DEVELOPMENT 

The  computer  program  was  written  in  Fortran  IV  com- 
puter language  using  double-precision  arithmetic  on  an 
IBM  360/67  digital  computer.   It  consists  of  a  main  calling 
program  and  seven  subroutines.   The  main  program  calls  the 
subroutines,  where  the  actual  calculations  occur,  in  the 
necessary  order.   Information  is  passed  between  the  main 
program  and  the  subroutines  by  the  use  of  a  common  storage 
core.   Appendix  E  contains  the  program  listing. 

A.   SYSTEM  DESCRIPTION 

After  the  piping  system  has  been  modeled  and  subdi- 
vided, subroutine  INPUT  is  used  to  read  the  data  input  of 
the  geometry  and  material  properties  of  the  system.   The 
geometry  of  the  piping  system  model  is  determined  by  the 
global  coordinates  of  the  nodes  of  the  model.   The  global 
system  is  a  right-handed  orthogonal  coordinate  system. 
Each  pipe  segment  must  be  parallel  to  a  global  axis,  but 
there  are  no  restrictions  on  the  placement  of  the  origin  of 
the  global  system.   The  parameters  from  the  subdivision  of 
the  system  are  read.   For  each  node  the  global  coordinates 
and  boundary  conditions  are  read.   If  the  node  is  fixed, 
the  stiffness  and  inertial  matrix  components  due  to  the 
node  displacements  are  not  assembled  into  the  system 
matrices . 
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The  numbering  of  nodes  and  elements  is  arbitrary. 
There  are  two  nodes  per  element  and  six  degrees  of  freedom 
per  node.   For  each  element,  the  node  numbers  and  pipe 
group  are  read.   A  pipe  group  includes  all  elements  that 
have  the  same  Young's  modulus,  Poisson's  ratio,  specific 
weight,  radius  of  gyration,  and  cross-sectional  dimensions. 
Appendix  D  shows  the  method  used  to  calculate  the  effect 
of  lagging  on  the  specific  weight  and  the  radius  of 
gyration. 

B.  STIFFNESS  AND  INERTIAL  MATRIX  FORMULATION 

Subroutine  F0RM  generates  the  element  stiffness  and 
inertia  matrices  as  shown  in  Appendix  A.   Rotatory  inertia 
and  shear  deformation  due  to  bending  are  neglected.   The 
orientation  of  the  element  is  determined  with  respect  to 
the  global  coordinate  system  and  the  correspondence  between 
local  and  global  degrees  of  freedom  is  established.   After 
determining  the  correspondence  between  the  local  and  global 
degrees  of  freedom,  the  element  stiffness  and  inertia 
matrices  are  assembled  to  form  the  system  stiffness  and 
inertia  matrices.   If  pipe  hangers  are  present,  the  corres- 
ponding (uniaxial)  hanger  stiffnesses  are  added  to  the 
appropriate  diagonal  elements  of  the  system  stiffness 
matrix. 

C.  FREE  VIBRATION  MODE  SHAPES  AND  FREQUENCIES 

Subroutine  CHJ0M0D  uses  Cholesky  decomposition  of  the 
system  inertia  matrix  and  coordinate  transformation  to  form 
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a  single  symmetric  matrix  for  which  eigenvalues  and  eigen- 
vectors are  found.   Subroutine  JACR0T  is  a  modification  of 
a  program  originally  coded  by  Professor  G.  Cantin.   It 
uses  the  Jacobi  variable  threshold  method  to  find  eigen- 
values and  eigenvectors  of  a  real  symmetric  matrix.   The 
modal  matrix,  when  transformed  back  to  system  coordinates, 
is  normalized  with  respect  to  the  system  inertia  matrix. 

Since  the  system  matrices  increase  in  size  with  finer 
subdivisions  or  more  complex  systems,  an  economizing  tech- 
nique, to  reduce  the  number  of  degrees  of  freedom  in  the 
eigenvalue  problem,  was  investigated  during  the  program 
development.   The  component  mode  synthesis  method  as  used 
by  Benfield  and  Hruda  [11]  was  studied.   It  was  established 
that  the  constrained  component  branch  technique  [11],  with 
interface  loading,  with  no  node  suppression  gave  results 
identical  with  those  obtained  by  the  direct  method  described 
above  when  applied  to  planar  vibration  of  a  clamped-c lamped 
beam  represented  by  six  elements.   Despite  the  success  of 
this  trial,  it  was  concluded  that  the  added  program  com- 
plexity accompanying  the  use  of  component  mode  synthesis 
would  offset  the  potential  economies.   Another  standard 
economizer  technique  [14]  was  considered,  but  ultimately 
rejected  for  the  same  reason. 

D.   MODAL  RESPONSE 

The  number  of  modes  to  be  used  in  the  stress  analysis 
and  the  specification  of  the  shock  spectrum  are  input  data 
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for  subroutine  SPCTRM.   There  are  three  choices  available 
to  select  the  number  of  modes  to  be  used:  a  designated 
upper  frequency  limit,  an  upper  limit  which  is  a  multiple 
of  the  fundamental  frequency,  and  an  option  for  any  method 
the  user  desires.   Three  choices  are  also  available  for 
specifying  the  shock  spectrum:  a  spectrum  defined  by 
straight-line  segments;  a  constant  velocity,  then  constant 
acceleration  shock  spectrum;  and  an  option  for  any  method 
the  user  desires.   The  shock  input  is  restricted  to  a 
translation  of  the  base  along  a  global  axis  and  is  speci- 
fied by  one  shock  spectrum  and  three  scaling  factors  for 
the  three  input  directions.   Subroutine  M0DE  modifies  the 
spectrum  velocity  for  modal  mass,  if  desired.   There  are 
three  choices  available:  no  correction,  log  of  the  modal 
weight  correction,  and  an  option  for  any  method  the  user 
desires.   This  subroutine  also  calculates  the  mode 
participation  factors. 

E.   STRESSES 

Subroutine  STRESS  calculates  stresses  at  the  nodes  of 
the  elements.   The  element  peak  displacements  are  deter- 
mined and  realigned  from  global  degrees  of  freedom  back  to 
the  local  element  degrees  of  freedom.   The  element  gener- 
alized force  vector  is  then  calculated.   Consider  the 
action  of  the  generalized  force  vector  on  the  element 
shown  in  Pig.  3- 
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FIG.  3 
STRESS  LOCATION  POINTS 

At  locations  a,b,c,  and  d,  the  stress  vector  acting  on  the 
cross-section  may  be  resolved  into  normal  and  shearing 
stresses.   Fig.  k   shows  the  positive  sign  convention  used 
for  the  normal  and  shearing  stresses. 

The  normal  stress  a  and  the  shearing  stress  x  at  each 
point  shown  in  Fig.  3  are  given  by: 
1.   at  location  a) 


a  =  -f-j/A  -  fV  D/2I 


x  =  -fjj  D/4I  +  f2  Q/2It 
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FIG.  4 
POSITIVE  STRESS  CONVENTION 

where  I  is  the  second  moment  of  cross-sectional  area  about  a 
diameter,  D  is  the  outside  diameter,  Q  is  the  first  mo- 
ment, about  a  diameter,  of  the  portion  of  the  cross-sectional 
area  lying  on  one  side  of  the  diameter,  and  t  is  the  pipe 
wall  thickness. 

2.   at  location  c) 


a  =  f  /A  +  fn  D/2I 


x  =  f1Q  D/4I  -  fg  Q/2It 


3.      at    location  b) 


a    =    -    f1/A   -    fg    D/2I 


f ..    D/JJI    +    f.    Q/2It 
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at  location  d) 


a   =  -  f?/A  -  f12  D/2I 


t  =  f1()  D/4I  -  f9  Q/2It 


These  normal  and  shearing  stresses  are  then  used  to  find 
the  octahedral  shearing  stresses  [12],  which  are  given 
by 


Toct  =  1/3  (2*2  +  6^)2  W 


For  each  element  the  above  calculations  are  performed  for 
each  mode.   At  each  of  the  four  points,  the  shearing 
stresses  for  the  individual  modes  are  then  combined  by 
determining  the  square  root  of  the  sum  of  their  squares. 
This  process  is  accomplished  for  each  input  direction. 
This  is  done  for  each  element  in  turn,  starting  with  the 
first. 

F.   PROGRAM  OUTPUT 

All  input  information  is  echo-checked.   The  square 
of  the  mode  frequency  and  the  mode  shape  are  printed  for 
each  mode.   The  spectrum  velocities  and  mode  participation 
factors  are  also  printed.   The  octahedral  shearing  stresses 
at  points  a,  b,  c,  and  d  (Fig.  3)  of  each  element  are 
printed  for  the  three  shock  input  directions. 
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IV.   PROGRAM  TESTING 

In  order  to  explore  the  capabilities  of  the  program 
and  verify  its  integrity,  a  number  of  plane  test  configur- 
ations were  studied.   These  are  described  in  the  following 
sections. 

A.   TEST  PROBLEM  1 

The  configuration  shown  in  Fig.  5  has  three  identical 
uniform  runs  of  equal  length  between  the  central  junction 
and  the  clamped  edges. 


/& 
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FIG.  5 
TEST  PIPING  SYSTEM  1 

The  mode  shapes  and  frequencies  for  in-plane  vibration  of 
this  system  were  studied  extensively  with  a  developmental 
program.   The  present  program  gave  identical  results  for  the 
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in-plane  modes.   The  out-of-plane  modes  exhibited  the 
required  symmetric  or  anti-symmetric  form. 

Additional  program  tests  on  this  configuration  involved 
different  element  numbering  and  node  numbering  and  six 
different  orientations  of  the  global  axes  relative  to  the 
structure.   Mode  shapes  and  frequencies  were  unaffected  by 
these  changes. 

B.   TEST  PROBLEM  2 

This  system,  shown  in  Fig.  6,  consisted  of  a  straight 
uniform  pipe  clamped  at  the  ends  and  represented  by  four 
elements. 
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FIG.  6 

TEST  PIPING  SYSTEM  2 

The  frequencies  found  for  the  lower  modes  showed  good 
agreement  with  exact  results.   Table  I  shows  the  comparison 
of  the  first  three  bending  frequencies  of  the  finite  element 
solution  with  exact  results  [13].   All  the  eigenvectors 
showed  the  required  symmetric  or  anti-symmetric  form. 
Bending  modes  occurred  in  pairs  having  equal  eigenvalues 
and  with  the  corresponding  eigenvectors  representing 
deflections  in  orthogonal  planes. 


31 


Stresses  found  for  this  system  were  studied  in  detail. 
At  each  node,  the  stresses  are  calculated  at  two  points  of 
the  cross-section.   Except  at  the  clamped  ends,  there  are 
two  adjacent  elements  sharing  each  node,  so  that  stresses 
at  these  sections  are  calculated  twice.   Complete  agreement 
was  found  between  these  two  sets  of  stress  evaluations. 
Also,  identical  stresses  resulted  from  shocks  in  the  two 
directions  perpendicular  to  the  length. 

TABLE  I 

COMPARISON  OF  FIRST  THREE  BENDING  FREQUENCIES  OF  TEST 
PIPING  SYSTEM  2  WITH  EXACT  RESULTS 


ode 

Exact 
(Hz) 

Finite 

Element 

(Hz) 

Percent 
Difference 

1 

422.32 

422.65 

0.08 

2 

1164.15 

1174.26 

0.9 

3 

2282.20 

2329.64 

2.1 

C.   TEST  PROBLEM  3 

The  equal- legged  configuration  of  Fig.  7  consists  of 
uniform  pipe  throughout.   Mode  shapes  of  this  system  again 
showed  the  required  symmetric  or  anti-symmetric  form. 
Likewise,  the  shock-induced  stresses  exhibited  the  expected 
symmetry  properties. 
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FIG.    7 
TEST   PIPING   SYSTEM   3 
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V.   EXAMPLE  PROBLEM 

A.  PIPING  SYSTEM 

The  example  system  analyzed  is  a  model  of  the  main 
steam  piping  system  on  a  modern  naval  vessel.   The  main 
steam  piping  from  the  boiler  to  the  rigid  cross-connect 
anchor  assembly  is  modeled.   The  ends  of  the  piping  system 
are  clamped,  and  there  are  two  hangers.   The  system  is 
three-dimensional  with  no  branches.   The  pipe  properties 
are: 

outside  diameter      7.625  inches 

inside  diameter       6.011  inches 

Poisson's  ratio       0.300 

Young's  modulus       30x10  psi 

specific  weight       0.327  lb/in3 

radius  of  gyration    4 . 04  inches 

The  specific  weight  is  a  fictitious  value  which  accounts 
for  the  weight  of  five  inches  of  lagging  bonded  to  the  pipe 
Fig.  8  is  a  schematic  of  the  piping  system.   The  piping 
section  PI  is  84  inches  long,  P2  is  264  inches  long,  P3  is 
108  inches  long,  and  P4  is  288  inches  long. 

B.  DATA  INPUT  AND  OUTPUT 

The  input  data  cards  are  punched  in  accordance  with  the 
instructions  contained  in  Appendix  E.   The  echo-check  of 
the  data  is  shown  in  Table  II  and  is  arranged  as  it  would 
appear  at  the  end  of  the  program  deck.   The  system 
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FIG.  8 
EXAMPLE  PIPING  SYSTEM 

eigenvalues  and  eigenvectors  are  listed  in  Appendix  C, 
Table  VIII.   Tables  III  and  IV  show  the  modal  velocities 
and  the  mode  participation  factors.   The  stress  output  is 
given  in  Table  V. 

This  problem  required  a  storage  capacity  of  17^K  bytes, 
and  2.67  minutes  of  computer  time. 
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TABLE   II      DATA   INPUT,  EXAMPLE   PROBLEM 
PROBLEM   NUMBER!  1 

NUMBER   OF    PROBLEMS    TO   BE    SOLVED 


TOTAL 
ELEMENTS 


12 


TOTAL 

NODES 


DEG   OF    FREEDOM 
PER    NODE 


DEG    OF    FREEDOM 
PER    ELEMENT 


12 


TOTAL 

GROUPS 


SYSTEM    TOTAL 
DEG    OF    FREEDOM 


66 


NODE 

X    COORD 

Y   COORD 

Z    COORD 

NODAL    BOUNDARY 

INCHES 

INCHES 

"  INCHES- 

CONDITIONS 

1 

0.0 

0 

0 

0.0 

J 

2 

42.00 

0 

.0 

0.0 

3 

84.00 

0 

0 

0.0 

0 

4 

84.00 

0 

0 

66.00 

0 

5 

84.00 

0 

.0 

132.00 

0 

6 

84.00 

0 

3 

198.03 

0 

7 

84.00 

0 

0 

264.00 

0 

8 

138.00 

0 

0 

264.00 

0 

9 

192.00 

0 

0 

264.00 

0 

10 

192.00 

72 

.00 

264.00 

0 

11 

192.00 

144 

00 

264.00 

0 

n 

192.00 

216 

00 

264.00 

0 

192.00 

288 

.00 

264.00 

ELEMENT 

LEFT 

RIGHT- 

GROUP 

NUMBER 

NODE 

NODE 

NUMBER 

1 

1 

2 

I 

2 

2 

3 

. 

3 

3 

4 

. 

4 

4 

5 

' 

. 

5 

5 

6 

6 

6 

7 

7 

7 

6 

8 

8 

9 

;. 

9 

9 

\l 

. 

10 

10 

11 

12 

12 

13 

OUTSIDE 

INSIDE 

POISSONS       SPECIFIC 

YOUNGS 

DIAMETER 

DIAMETER 

RATI 

MODULUS 

IN. 

IN. 

LB/CU.IN 

PS  I 

SHEAR 

MODULUS 

PSI 


7.625 


6.Q11 


0.300 


3.000D    07 


1.1540    07 


PIPE 

GROUP 


RADIUS    OF    GYRATION 
INCHES 


1  4.  04 

NUMBER    OF    HANGERS: 


GLOBAL    DOF 

IADOF 

HANGER    STIFFNESS 
LB/IN. 

7 
11 

3 
3 

9.817477D    05 
9.8174770    05 

SPECTRUM 
TYPE 

1 

FREQUENCY 
TYPE 

BREAK    FREQUENCY 
HERTZ 

CONSTANT    VELOCITY 
INCHES/SEC. 

50.00 

100.00 

SHOCK    INPUT 
X 

FACTORS 
Y 

Z 

1.000 

1.000 

1.000 

MODE   WEIGHT    CORRECTION    TYPE 
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TABLE  III   MODE 

VELOCITIES  EXAMPLE,  PR( 

MODE 

SPECTRUM  VELOCITY 

INCHES/SEC. 

1 

100.00 

2 

100.00 

3 

10  3.00 

4 

100.00 

5 

100.00 

6 

10  0.00 

7 

100.00 

8 

100.00 

9 

100.00 

10 

88.87 

11 

84.49 

12 

75.10 

13 

71.39 

14 

58.88 

15 

53.46 

16 

47.11 

17 

42.75 

18 

42.19 

19 

32.52 

20 

30.30 

21 

29.92 

22 

29.23 

23 

27.07 

24 

26.08 

25 

25.76 

26 

25.43 

27 

22.07 

28 

21.76 

29 

17.29 

30 

17.04 

31 

16.10 

32 

15.18 

33 

14.34 

34 

14.14 

35 

13.84 

36 

12.24 

37 

12.16 

38 

12.08 

39 

10.67 

40 

10.29 

41 

10.03 

42 

9.65 

43 

9.32 

44 

8.85 

45 

8.35 

46 

8.31 

47 

8.02 

48 

7.82 

49 

7.53 

5  3 

7.14 

51 

6.91 

52 

6.61 

53 

6.51 

54 

5.94 

55 

5.80 

56 

5.48 

57 

5.27 

58 

4.74 

59 

4.58 

60 

4.47 

61 

4.30 

62 

3.91 

63 

3.68 

64 

3.43 

65 

2.97 

66 

2.73 
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TABLE  IV  MODE  PARTICIPATION  FACTORS,  EXAMPLE  PROBLEM 

MODE  MODE  PARTICIPATION  FACTORS 

X  Y  Z 


1 

2.42DD     DO 

-5.433D-01 

-2.376D-01 

2 

8.881D-01 

1.936D    00 

-8.579D-03 

3 

6.749D-01 

-2.151D-01 

1.366D    00 

4 

1.983D-01 

-8.6390-01 

-2.739D-01 

5 

1.813D-01 

3.681D-02 

4.449D-01 

6 

1.0D2D     JJ 

2.950D-02 

-4.5250-01 

7 

-2.002D-01 

8.847D-01 

-5.983D-02 

8 

1.342D-01 

2  .9090-02 

-4.1320-01 

9 

3.281D-01 

7.394D-02 

-1.452D-01 

10 

-5.372D-01 

2.087D-02 

-2.2340-01 

11 

-1.760D-01 

-2. 7650-02 

3.701D-01 

12 

-1.461C- 01 

1.7600-02 

-4.3580-01 

13 

-2.436D-02 

-2.057D-02 

-2.5810    00 

14 

2.101D-01 

6.4130-01 

-2.771D-02 

15 

2.120D-01 

3.6170-02 

6.2730-02 

16 

2.414D-01 

-4.2060-01 

-5.3080-02 

17 

-1.557D-01 

-1.836D-01 

-8.9290-03 

18 

-2.999D-01 

2.7700-01 

-9.8480-02 

19 

-2.523D-03 

-1  .5750    00 

2.7830-02 

20 

3.7280-02 

3.9690-01 

4.6230-01 

21 

-7.001D-02 

-7.5210-01 

1.8320-01 

22 

•      5.6800-02 

-2.8310-03 

3.3900-01 

23 

2.730D-02 

3.0590-02 

1.122D-01 

24 

1.372D-02 

-1.4560-01 

-7.496D-02 

25 

2.O96D-01 

3.593D-02 

5.776D-02 

26 

-5.284D-02 

1.7100-01 

-2.951D-03 

27 

3.018D-01 

1.8490-02 

-2.0460-01 

28 

6.857D-02 

-4.139D-01 

3.6120-32 

29 

1.502D-01 

1.9650-02 

2.480D-01 

30 

1.343D-01 

-3.7780-02 

5.7400-02 

31 

1.741D-01 

-5.5900-02 

-1.286D-01 

32 

8. 7360-02 

-7.8550-03 

-2.3340-01 

33 

3.208D-01 

-3.3680-02 

-1.143D-01 

34 

2.O57D-01 

1.0940-01 

-7.9530-02 

35 

-1.832D-02 

1.4720-03 

-1.7040-01 

36 

-5.653C-02 

-2.6610-01 

3.9710-02 

37 

-2.768D-01. 

-2.8490-02 

-1.585D-01 

38 

-1.338D-01 

4.8110-02 

1.0380-01 

39 

5.2820-31 

1.843D-02 

-1.1140-01 

40 

-3.848D-01 

-1.9360-01 

-4.5490-02 

41 

3.217D-01 

-4.5550-01 

1.7640-03 

42 

-9.8130-02 

-2.3620-01 

5.5900-02 

43 

1.810D-01 

9  .5520-02 

8.5750-02 

44 

-4.923D-31 

-6.4730-03 

-6.5880-02 

45 

-1. 517D-01 

4.6710-03 

-9.275D-02 

46 

4.5190-02 

2.3120-01 

6.8950-C3 

47 

-3.715D-32 

8.8880-04 

-4. 6700-02 

48 

-7.110D-02 

2.7250-02 

-1.0660-01 

49 

-8.744D-02 

4.0900-02 

2.3070-02 

53 

-5.489D-02 

-6.571D-02 

-1. 5690-02 

51 

3.389D-02 

4.8980-02 

-3.5650-02 

52 

8.851D-02 

-1.2770-01 

5.838D-02 

53 

-5.329D-02 

-6.4960-02 

-1.231D-01 

54 

5.C54D-02 

1  .6870-01 

8.028D-03 

55 

6.4810-34 

-2.6370-02 

-1.6450-02 

56 

-2.0470-02 

2.8140-01 

-7.3650-03 

57 

-2.748D-03 

-5.7720-02 

-3.583D-03 

58 

1.960D-02 

-2.5910-02 

-1.7730-02 

59 

-1.5410-C1 

1.049^-02 

8.5290-02 

60 

-1.5480-01 

-1.8900-02 

1. 1470-01 

61 

-8.0590-03 

2.6220-02 

4.9190-02 

62 

1.948D-03 

-1.546D-01 

2.3570-03 

63 

8. 5170-05 

-1.3750-01 

1  .0660-04 

64 

-3.405D-02 

-5.9190-03 

-3.0  030-02 

65 

5. 4040-03 

-7.4280-03 

6.0630-03 

66 

2.  524D-01 

4. 6170-35 

3.8080-02 
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TABLE  V     OCTAHEDRAL   SHEARING 
ELEMENT  1 

9.271550    03      6.508980    03      8, 

6.674260   03       1.687530    04      6. 

6.609590    03      6.198620    33      2. 

4.588570   03      8.645530    03       3. 


STRESSES,  EXAMPLE   PROBLEM,    PSI 

jl    ELEMENT  7 

408350   03  ^        1.055150   04      3.19562D    03  7.845000    03 

48256D    03^       4.203340    03       5.39532D    03  2.892360    03 

61721D    03  5.88687D   03       3.140880    03  3.869280    03 

095510    03         4.372620    03      5.55352D    03  3.82004D    03 


ELEMENT  2 

6.609590   03  6.198620   03  2. 

4.58857D    03  8.64553D    03  3. 

1.879580   04  7.527770    03  7. 

4.96003D   03  7.392360    03  1, 


j  ELEMENT           8 

617210    03 j  5.886870    03 

095510    03  4.37262D   03 

71406D   03  2.656300    03 

04808D    03  7.424390    03 


3.14088D  03 

5.55052D  03 

3.213450  03 

1.318900  04 


3.869280  03 
3.82004D  03 
3.567750  03 
4.40520D   03 


ELEMENT  3 

5.764590  03 
1.88813D  04 
4.458580  03 
1.114120    04 


7.613430  03 
6.471280  03 
5.76172D  03 
4.91469D    03 


1.77620D  03 
7.52389D  03 
3.03388D  03 
2.67322D    03 


ELEMENT  9 

7.59211D  03 
2.46920D  03 
5.77644D  33 
4.85717D   03 


1.27242D    04  4.36537D   03 

2.81124D    03  3.50809D    03 

1.01342D    04  3.201110    03 

2.06310D    03  8.55130D    03 


ELEMENT  4 

4.458580  03 
1.11412D  04 
4.561880  03 
7.49948D    03 


5.761720  33 
4.91468D  03 
9.146560  03 
4. 300960    33 


3.03388D  03 
2.673220  03 
3.568610  03 
6.18924D    03 


ELEMENT  10 
5.77644D  03 
4.857170  03 
7.21896D  03 
9.86817D    03 


1.01042D   04  3.20111D    03 

2.06310D    03  8.55130D    03 

6.12863D    03  4.520530    03 

3.488530    03  1.570740    04 


ELEMENT  5 

4.561880  03 
7.499480  03 
4.543860  03 
6.25280D   03 


9.146560  03 
.4.300960  03 
7.59655D  03 
4.33059D    03 


3.56861D  03 
6.18924D  03 
4.76524D  03 
7.081600    03 


ELEMENT  11 
6.94094D  03 
9.868170  03 
1.04300D  04 
3.535350    03 


6.083360    03  4.586790    03 

3.488530    33  1.573740   04 

4.37922D    03  2.33575D    03 

2.30537D    03  9.55394D    03 


ELEMENT  6 

4.543860  03 
6.252800  03 
3.61390D  03 
1.086730   04 


7.59656D  03 
4.33059D  03 
4.921470  03 
4.381410    03 


4.765240  03 
7.08160D  03 
2.47225D  03 
7.711130  03 


ELEMENT  12 

1.34303D  04 

3.535340  03 

1.769770  04 

6.153500  03 


4.379220    33  2.33575D    03 

2.30537D    03  9.553940   03 

1.275160    04  6.00477D    03 

2.751500    03  1.630620   04 
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VI.   DISCUSSION 

Some  of  the  limitations  of  the  capabilities  of  the 
program  developed  above  are  considered  here. 

A.  CONFIGURATION  LIMITATIONS 

The  restrictions  that  pipe  axes  must  be  parallel  to  a 
global  axis  and  that  finite  radius  bends  are  not  represented 
provide  the  principal  limitations  on  the  configurations  that 
can  be  modeled.   There  are  no  inherent  limitations  on  allow- 
able topological  complexity. 

Certain  additional  features  of  real  piping  systems  that 
are  not  represented  in  the  present  modeling  include  added 
mass  due  to  fittings  and  pipe  contents,  and  partial  fixity 
or  elastic  restraint  at  ends  or  intermediate  points. 

Despite  these  limitations,  it  is  believed  that  a  signi- 
ficant fraction  of  current  piping  systems  can  be  adequately 
modeled  using  the  present  program. 

B.  SIZE  LIMITATIONS 

For  the  present  purpose,  the  appropriate  measure  of 
system  size  is  the  number  of  degrees  of  freedom  of  the 
system  model.   This  number  n,  which  bears  no  direct  relation 
to  physical  size,  determines  the  core  storage  requirements 

and  execution  time  of  the  program.   For  approximate  estima- 

2 
tion,  storage  requirements  are  proportional  to  n  and 

execution  time  is  proportional  to  n  .   Using  the  data  from 

the  example  problem,  one  can  estimate  that  a  100 
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degree-of-freedom  system  would  require  about  400K  bytes  of 
core  storage  and  about  9  minutes  execution  time. 

Increasing  the  problem  size  also  increases  the  round-off 
errors.   Because  double-precision  arithmetic  (56  bit  mantissa) 
is  used  throughout,  observed  round-off  effects  have  been 
found  negligible  (a  maximum  of  1  unit  in  the  6th  significant 
digit  for  n  ■  66)  in  all  applications  to  date.   The  very 
small  errors  detected  are  attributed  to  the  residual  eigen- 
vector errors  present  when  the  Jacobi  rotations  are 
terminated. 

In  view  of  the  foregoing  considerations,  it  is  believed 
that  the  practical  upper  limit  on  problem  size  (with  the 
IBM  360/67)  is  about  n  =  110  and  is  determined  by  core 
storage  capacity. 


Ml 


VII.   CONCLUSIONS 

It  is  concluded  that  an  effective  program  has  been 
developed  for  determining  shock-induced  stresses  in  piping 
systems.   The  principal  limitations  of  the  present  version 
can  be  removed  by  adding  features  that  are  clearly  within 
the  current  state-of-the-art.   Recommended  additions  are 
listed  below: 

1.  Remove  the  pipe  axes  orientation  restriction  so  that 
a  general  piping  system  can  be  analyzed. 

2.  Develop, element  stiffness  and  inertia  matrices  for 
bends. 

3.  Modify  the  program  so  that  partial  fixity  and  elastic 
restraint  at  the  ends  or  intermediate  points  may  be 
included  in  the  boundary  conditions. 

4.  Include  the  mass  effects  of  valves,  fittings,  and 
pipe  contents. 

5.  Replace  the  Jacobi  rotation  method  by  a  more  efficient 
eigenvalue-eigenvector  algorithm. 
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APPENDIX  A 
ELEMENT  STIFFNESS  AND  INERTIA  MATRICES 
TABLE  VI  Element  Stiffness  Matrix 
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APPENDIX  B:   VELOCITY  SHOCK  SPECTRA  AND  MODAL  RESPONSE1 


A  shock  spectrum  exhibits  the  maximum  response  dis- 
placement of  a  single  degree-of-freedom  system  whose  base 
is  subjected  to  the  shock  motion.   Consider  a  single  degree- 
of-freedom  whose  base  has  a  shock  displacement  s,  and  whose 
displacement  relative  to  the  base  is  z.   The  equation  of 
motion  is 


z  +.  0)  z  =  -  s  (B.l) 


where  u  is  the  natural  circular  frequency  of  free  vibration 

with  the  base  fixed  (s  =  0).   If  z    represents  the  ex- 

max 

treme  value  of  z  in  response  to  the  shock  motion  s,  then  the 

displacement  shock  spectrum  for  that  shock  motion  is  a  plot 

of  z    versus  w.   Equivalent  information  can  be  given  in 
max 

the  form  of  a  velocity  shock  spectrum.   In  this  form,  the 
spectrum  velocity  V  is  related  to  the  maximum  displacement 

zmax  by 


max 


To  develop  the  equations  for  finding  modal  response 
from  shock  spectrum  data,  let  w  be  the  vector  of  absolute 


This  Appendix  is  based  on   material  presented  by  Prof 
R.E.  Newton  in  the  course  ME4522,  September  1971  - 
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displacements  in  the  system  degrees-of-freedom.   For  a 
uniaxial  translation  s  of  the  base,  this  may  be  expressed 
as 


C[  +  us 


(B.3) 


where  q  is  the  vector  of  displacements  relative  to  the 
base,  and  u  is  a  vector  of  constants.   To  illustrate  the 
meanings  of  these  vectors,  consider  the  beam  element  of 
Fig.  9. 


\\\   Inerital  Reference        \\  \\ 


FIG.  9 
BEAM  ELEMENT  COORDINATES 


m  m 

For  this  example,  q  =  [qn  q2  q^  q^]  ,  w  =  [q1+s  q2  q^  +  s  q^]  , 
and  u  =  [1  0  1  0].   It  can  be  seen  that  u  represents  the 
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absolute  displacements  resulting  from  unit  base 
translation. 

The  equations  of  motion  of  an  n  degree-of-freedom 
system  with  base  motion  may  be  written 

M  w  +  K  q_  =  0 

Using  Eq.  B.3,  this  becomes 

Mqt+K^=  -Mus 

Using  the   substitution  q  =  V  p_  where   p   is   the   vector   of 
principal   coordinates,    this  may  be   rewritten  as 


VTMV£+VTKVp=   -VT  Mus 


or 


p_+fi2p   =   -bs  (B.4) 


T 
where  b  =  V  M  u.   Eq.  B.4  is  equivalent  to  the  scalar 

equations 


Pr  +  "r2  Pr  =  "  br  s  ,   (r=l,2,...,n)       (B.4a 


where  p   and  b   are  the  rth  components  of  p_  and  b, 
respectively,  and  w   is  the  circular  frequency  of  mode  r 
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The  component  b   Is  called  the  mode  participation  factor 
for  mode  r.   It  is  given  by 


T 
br  =  v(r)   M  u  (B.5) 


where  v    is  the  eigenvector  for  mode  r  (rth  column  of 

I)- 

Comparing  Eq.  B.4a  with  Eq.  B.l,  and  using  Eq.  B.2, 

the  maximum  response  in  mode  r  may  be  expressed  as 


(p  )     =  b   V  /oi  (B.6) 

Vir;max    r  r  r 


where  V  is  the  spectrum  velocity  at  frequency  u  .   Using 


r 


Eq.  B.6,  the  corresponding  extremom  for  the  relative 
displacement  vector  q_  is 


q    =  v(r)   b   V  /u)  (B.6a) 

^max   -      r  r  r 


Eq.  B.6a  is  used  for  calculating  peak  modal  responses  in 
order  to  find  stresses. 

In  shipboard  shock  studies,  it  has  been  found  that 
massive  equipment  may  partially  suppress  the  base  motion 
[6].   This  may  be  taken  into  account  by  appropriate  modifi- 
cation of  the  spectrum  velocity  V  .   Such  a  correction  re- 
quires apportionment  of  the  system  mass  among  the  modes. 
The  modal  mass  for  the  rth  mode  is  taken  to  be 


m  =  b 
r    r 


2  (B.7) 
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APPENDIX  C:  EXAMPLE  PROBLEM  EIGENVALUES  AND  EIGENVECTORS 

Table  VIII  shows  the  first  six  eigenvalues  and  the 
corresponding  eigenvectors  of  the  example  problem.   The 
eigenvalue  (square  of  the  circular  frequency)  appears  at 
the  head  of  a  column  followed  by  the  components  of  the 
eigenvector  grouped  by  nodes.   The  remaining  eigenvalues 
and  eigenvectors  are  omitted.   Spacing  of  the  eigenvalues 
remains  approximately  uniform  throughout  the  remainder  of 

o 

the  set.   The  largest  eigenvalue  is  1.32  x  10  . 
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TABLE  VIII 
EXAMPLE  PROBLEM  EIGENVALUES  AND  EIGENVECTORS 
3.20210D  02    1.768220  03    4.283940  03    7.284420  03 


1.366430  04    1.93694D  04 


3.734190-05 
2.41956D-03 
4.49829D-03 
2.61317D-04 
-9.722550-05 
1. 227360- 04 


1.34663D-05 
1.09003D-02 
1.04163D-03 
•1.29764D-03 
■2.247990-05 
3.964450-04 


2. 084720-04 
■2.65694D-02 
1.  106640-02 
4.11111D-04 
-2.55767D-04 
■1.134310-03 


-1.14030D^04 
-1.28407D-01 
-5.  213800-03 
1.33441D-03 
1  .123570-04 
-5.43247D-03 


■4.651040-04 
8.53335D-02 

•1  .59767D-02 
1.670380-04 
3.099150-04 
3.61351D-03 


9.202400-04 
5.243660-02 
2.80391D-02 
3.87608D-04 
•5.48334D-04 
2.206790-03 


7.  4  6  83  00-  05 

1.09429D-02 
■1.65671D-03 
5.226180-04 
5. 072280-04 
2.90708D-04 


2.693140-05 
2.303660-02 
■3.  865930-04 
■2.59488D-03 
1.175200-04 
6.000790-05 


4. 16900D-04 
■8.44930D-02 
■1.211330-03 
8. 219070-04 
1.  109210-03 
■1.50032D-03 


-2.28018D-04 
■4.007280-01 
1.910950-03 
2.66708D-03 
■5.85658D-04 
-6.91251D-03 


-9.29892D-04 
2.67855D-01 
1.14491D-02 
3.336680-04 

■2.05621D-03 
4.70215D-03 


1.83959D-03 
1.628400-01 
■1.906140-02 
7.738690-04 
3.551720-03 
2.82479D-03 


7.586450-02 
-3.27225D-02 
-1.502070-03 
7.81717D-04 
1.68048D-03 
6.80364D-04 


1.66544D-02 
2.268140-01 
-3.50837D-04 
•3.29915D-03 
3.47592D-04 
-1.36944D-03 


1. 461590-01 
■1.37883D-01 
-8.567620-04 
5.61288D-04 
2.712830-03 
-1.50438D-03 


■7.41273D-02 
-5. 514730-01 
1.733630-03 
9.10543D-04 
■1.33431D-03 
-6. 234990-03 


•2.396660-01 
2.05302D-01 
1.086020-02 
1.  877200-03 

-3.91045D-03 
4.76584D-03 


935970-01 
26126D-02 
804080-02 
71397D-03 
854290-03 
827700-03 


2.077980-01 
■8.98682D-02 
■1.34738D-03 
9.335880-04 
2.214390-03 
1.06992D-03 


4.113490-02 
4.263500-01 
■3.150070-04 
-2.55823D-03 
3.624590-04 
-2.79781D-03 


2.873100-01 
-1.35861D-01 
■5.01745D-04 
■7.25499D-04 

1. 113250-03 
•1.505580-03 


-1.40471D-01 
-4.479070-01 
1.55477D-03 
-4.14818D-03 
-4.794290-04 
■5.5373  5D-03 


■3.879600-01 
2.414480-02 
1.025300-02 
3.312280-03 

-4.39684D-05 
4.800720-03 


5.383140-01 
■4.12832D-02 
■1.697730-02 

1.745820-03 
■2.271960-03 

2.80638D-03 


3.55428D-01 
-1.54172D-01 
■1.  192630-03 
1.005280-03 
2.171520-03 
1.45933D-03 


6.127200-02 
5.49921D-01 
•2.791070-04 
■1.  16814D-03 
2.32548D-04 
-4.223990-03 


2.50525D-01 
•3.75462D-02 
■1.46464D-04 
•2.229600-03 
■2.408720-03 
■1.  503930-03 


■1.256530-01 
■4.37906D-02 

1.37451D-03 
•7. 483740-03 

8.999310-04 
■4.821860-03 


-2.50044D-01 
■1.78932D-01 
9.62857D-03 
2.470720-03 
3.612420-03 
4.806570-03 


1.421630-01 
■1.13515D-01 
■1.58733D-02 

3.151950-04 
-8. 321710-03 

2.761020-03 


4.  84437D-01 
-2.219640-01 
•1.037840-03 
1.050640-03 
1.67614D-03 
1.84853D-03 


7.  136440-02 
5. 922980-01 
-2.43147D-04 
-3.  187010-04 
8.26620D-05 
-5.64686D-03 


•3.727910-02 
1.513220-01 
2.088940-04 
•3.438590-03 
•6.26385D-03 
-1.49942D-03 


-3.99987D-02 
4.44113D-01 
1.19302D-03 

■6.623490-03 
1. 465560-03 

-4. 090820-03 


-5.137650-02 
-2.  818790-01 
8.98798D-03 
8.15115D-04 
1.02998D-03 
4.783360-03 


■3.036070-01 
-9.597510-02 
■1.473150-02 
-5.005620-04 
-3.163490-03 
2.69201D-03 


4.844510-01 
-1.159680-01 
-7.857740-02 
1.11158D-03 
1. 227640-03 
2.06739D-03 


-5 


7.13490D-02 
2.76427D-01 
2.70496D-03 
1.87342D-04 
1.566670-05 
781560-03 


■3.73090D-02 
7.10119D-02 
3.940520-01 

■4.58183D-03 
■7. 87601D-03 
-1.448660-03 


•3.98607D-02 
2.25377D-01 
7.170640-02 

•3.766930-03 
1. 203200-03 
4.12090D-03 


-5.06787D-02 
-5.032210-02 
7.37544D-02 
3.146250-06 
-2. 844290-03 
3.23151D-03 


•3.04055D-01 
2.482560-02 
■3.788760-02 
•7.11803D-06 
3.160610-03 
1.32777D-03 


4.844510-01 
■5.323830-05 
■1.372120-01 
1.  172420-03 
9.79064D-04 
2.211200-03 


7.13231D-02 
1.192160-03 
-2.735590-03 
-5.58841D-05 
-6.817080-06 
-4.079800-03 


•3. 732580-02 
2.462120-04 
8. 116120-01 
-5.71925D-03 
-7.355270-03 
•1. 10736D-03 


-3.969890-02 
4.44083D-04 

-1.282270-01 

-9.02245D-04 
9.235360-04 

-4.11453D-03 


•4.99239D-02 
■2.334730-03 
2.579220-01 
■8.08835D-04 
-3.532320-03 
-2.12272D-03 


-3.040190-01 
-1.372800-03 
-2.65755D-01 
4.863670-04 
4.543460-03 
-2.71595D-03 


3.223770-01 

-3.993090-05 
-5.555850-02 
1.023740-03 
7.344430-04 
2.254460-03 


2.42176D-01 
8.943870-04 
-2.962940-03 
3.44072D-05 
-5.118530-06 
-7.87964D-04 


1.77887D-02 
1.847940-04 
3.684660-01 
■6.240510-03 
•5.531030-03 
•4.405480-04 


2.072240-01 
3. 334770-04 

-1.034650-01 
1.293160-03 
6.95771D-04 

-2.32091D-03 


3.027290-01 
■1.75514D-03 

1.527010-01 
-2.121270-03 
■2.67172D-03 
-5.61151D-03 


6.34815D-02 
•1.03302D-03 
•1.749090-01 
2.211040-03 
3.44875D-03 
•6.110440-03 


1.67489D-01 
■2. 662160-05 
■6.52906D-04 
4.  239510-04 
4.89697D-04 
1.98080D-03 


2.105480-01 

5.963870-04 
-8.49562D-05 

2.89245D-05 
-3.415200-06 

1.42427D-03 


2.93284D-02 
1.232610-04 
5.64067D-03 
•3.012020-03 
•3.694300-03 
8.409130-05 


2.586930-01 
2.22516D-04 
-2.71326D-03 
9.958980-04 
4.653360-04 
8.950440-04 


5.251100-01 
-1.17205  0-03 

3.42056D-03 
-1.42807D-03 
■1.79190D-03 

3.088030-04 


3.673960-01 
-6.90308D-04 
•4.391680-03 
1.73252D-03 
2.31889D-03 
•1.28743D-03 


4.818110-02 
■1. 331100- 05 

7.315500-03 
-9.92107D-05 

2.448690-04 

1.24990D-03 


8.04306D-02 
2.98232D-04 
4.820690-04 
-6.354490-06 
-1.708410-06 
1.840420-03 


1.346290-02 

6.  164940-05 
•5.242910-02 

6.96480D-04 
■1. 84924D-03 

2.856850-04 


1.17881D-01 
1.113170-04 
1.714040-02 
-2.21904D-04 
2.331150-04 
2.488810-03 


2.76985D-01 
•5.866 110-04 
■2.559000-02 

3.24666D-04 
■8.991850-04 

5.464770-03 


2.298260-01 
•3.45644D-04 

3.17866D-02 
■3.929870-04 

1.16538D-03 

4.226120-03 
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APPENDIX  D:  CALCULATION  OF  MODIFIED  SPECIFIC  WEIGHT 
AND  RADIUS  OF  GYRATION 

There  is  negligible  effect  on  the  stiffness  of  the 
pipe  element  from  the  lagging  on  the  pipe;  however,  the 
lagging  does  contribute  significantly  to  the  mass  of  the 
element.   By  appropriately  modifying  the  pipe  element 
specific  weight  and  radius  of  gyration,  the  mass  effect  of 
the  lagging  can  be  included.   The  modified  specific  weight 
is 


y   =  (w-   +  WT  )/A  (C.l) 

m      p     j_i 


where  y      is  the  modified  specific  weight  of  the  pipe  ele- 
ment, and  the  sum  (W  +  WT )  is  the  combined  weight  of  the 
p    L 

lagging  and  pipe  per  unit  length.   The  radius  of  gyration 
is  given  by  the  relation 


J   +  JT  =  (ni  +  mT)  r2  (C.2) 

p     L      p     L 


where  J  and  JT  are  the  mass  moment  of  inertia  per  unit 
p      L 

length  of  the  pipe  and  lagging,  respectively,  (m  +  mT )  is 

p    L 

the  combined  mass  of  the  pipe  and  lagging  per  unit  length, 
and  r  is  the  radius  of  gyration. 
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APPENDIX  E:  PROGRAM  LISTING 
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